A Simple Model for Plastic Dynamics of a Disordered Flux Line Lattice 
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We use a coarse-grained model of superconducting vortices driven through a random pinning poten- 
tial to study the nonlinear current- voltage (IV) characteristics of flux flow in type II superconductors 
with pinning. In experiments, the IV relation measures flux flow down a flux density gradient. The 
work presented here treats this key feature explicitly. As the vortex repulsion weakens, the vor- 
tex pile maintains a globally steeper slope, corresponding to a larger critical current, for the same 
pinning potential. In addition, the magnitude of the peak in the differential resistance falls as the 
resistance peak shifts to higher currents. The model also exhibits so-called "IV fingerprints", and 
crossover to Ohmic (linear) behavior at high currents. Thus, many of the varieties of plastic be- 
havior observed experimentally for soft flux line systems in the "peak regime" are reproduced in 
numerical simulations of the zero temperature model. This model describes a two-dimensional slice 
of the flux line system at the scale of the London length (A). It does not include any degrees of 
freedom at scales much smaller than A, which are required to specify the degree of disorder in a flux 
line lattice. Instead, the nonlinear transport behaviors are related to the self-organized, large scale 
morphologies of the vortex river flow down the slope of the vortex pile. These morphologies include 
isolated filamentary channels, which can merge at higher flow rates to make a braided river, and 
eventually give uniform flow at even higher flow rates. The filamentary structure is associated with 
an IV characteristic that has concave, or positive, curvature. The braided river is associated with 
the peak in the differential resistance, where the curvature of the IV relation changes to convex. 
The transition to Ohmic behavior comes about as the braided river floods when it cannot absorb 
a higher level of flow. We propose that these self-organized morphologies of flux flow down a flux 
gradient govern the various plastic flow behaviors, including nonlinear IV characteristics, observed 
in type II superconductors with random pinning. 
PACS numbers: 74.60Ge, 74.60Jg, 64.60Ht, 62.20.Fe 



I. INTRODUCTION 

Collective transport in disordered media is a 
widespread and poorly understood phenomena. A great 
deal of experimental and theoretical effort in this area 
has been devoted to studying the nonlinear dynamics of 
the disordered flux line lattice (FLL) in type II super- 
conductors. The FLL exhibits a threshold behavior due 
to the competition between pinning and flux line repul- 
sion |jlg|. In response to a force, such as that associated 
with a transport current, the three dimensional FLL can 
move smoothly via elastic deformations, maintaining its 
integrity and order. However, in another regime the FLL 
deforms plasticly. In that regime disorder becomes more 
important, and the moving FLL manifold tears as some 
flux lines (in two-dimensions, vortices) move while oth- 
ers do not. As a result, the flow pattern breaks up in a 
nonuniform way. It is generally believed that the micro- 
scopic structure of the FLL, or defects in it, is fundamen- 
tal to transport behavior both in the plastic and elastic 
regimes. 

Part of the attention to the dynamics of a moving 
FLL has been motivated by interests in possible, exotic 
phase transitions and glassy phases, melting and other 



complicated scenarios associated with structural order 
in the FLL. Most experiments, however, are essentially 
transport studies, and are, quoting Higgins and Bhat- 
tacharya || , "notoriously ill-suited for the study of ther- 
modynamic phase transitions. These experiments yield 
direct information only about the mobility of flux lines, 
i.e. on dynamics and pinning, which would then have to 
be connected, through highly model-dependent ways, to 
the structure of the state they pertain to." 

Here, we show that many of the empirical results found 
in transport studies of plastic flux flow in type II super- 
conductors may be obtained with an extremely simple 
model B. It has been recognized for many years that, in 
the presence of pinning, magnetic flux in type II super- 
conductors forms a pile with an overall, global slope, akin 
to a sandpile. Penetration of magnetic flux into super- 
conductors driven solely by the flux density gradient has 
been described using molecular dynamics (MD) simula- 
tions [J5|— JT0|] , and by the model used here [0 . However, the 
flux gradient has not been taken into account in any pre- 
vious numerical simulation studies of the current-voltage 
characteristic. One possible reason is that previous nu- 
merical studies of flux motion at the scale of the vortex 
cores have not been able to reach a sufficiently large sys- 



tern size. The Bassler-Paczuski (BP) model j|, on the 
other hand, is a coarse-grained model and describes the 
magnetic flux dynamics at the much larger scale of the 
London length, making the large system size limit much 
more accessible. 

The numerical simulations presented here of the BP 
model 0] show that nonlinear behaviors, characteristic 
of experimental transport measurements of plastic flow 
in superconductors, arise as a result of vortex flow down 
a vortex density gradient. The effect of the transport cur- 
rent is modelled by a shift in boundary conditions, which 
leads to a generalized "tilt" of the vortex pile. Eventu- 
ally the "tilt" is sufficiently great so that some vortices 
can flow down the pile in the steady state, leading to the 
onset of a finite voltage. The vortex flow forms a vari- 
ety of river morphologies depending on the interaction 
strength between the vortices, compared to pinning, and 
the overall rate of flow. The flow patterns of magnetic 
flux are self-organized together with the magnetic flux 
profile, which is the substrate on which the flow takes 
place. Since the BP model does not contain any detailed 
information on the positions of the vortex cores at the 
micro-scale of the FLL, it cannot exhibit any structural 
ordering or disordering behavior. 

Although the BP model can also be studied at finite 
temperatures, here we use the zero temperature limit 
where thermal fluctuations of the flux motion may be 
ignored. The zero temperature approximation seems rea- 
sonable to describe the plastic transport dynamics of the 
low temperature superconductors and, perhaps, some as- 
pects of the high temperature superconductors as well. 
Also, we consider the limit where the depairing current 
density, jo, is extremely large compared to the critical 
current density, j c . This corresponds to the so called 
"weak pinning" regime. 

We argue that many of the varieties of collective 
transport dynamics observed in superconductors may be 
generic to repulsive particle systems driven through a 
disordered media. These behaviors are directly related 
in our coarse-grained model to large scale morpholo- 
gies of flow down a vortex density gradient and changes 
from filamentary strings, to a braided river, to uniform 
flow at high applied currents. Since the BP model ar- 
guably contains the essential physics of the disordered 
flux line system, at a coarse-grained level, and repro- 
duces a wide variety of experimental results on transport 
properties, we propose that these self-organized, large 
scale flow morphologies are also governing the nonlinear, 
plastic dynamics in the actual physical system: flux lines 
driven through a superconductor with a disordered pin- 
ning landscape by an applied transport current. 



A. Summary 

In the next section, we focus attention on an exten- 
sive study of 2H — NbSe2 as summarized by Higgins and 



Bhattacharya || . The reader may already compare Fig. 1 
and Fig. 2, which are schematic reproductions of experi- 
mental results, with Figs. 6 and 7, which are results from 
numerical simulations presented in this work. Section 
III presents a short review the standard theoretical pic- 
ture, associated with the work of Shi and Berlinsky pl| , 
on tearing of the FLL and the resultant nonlinear IV as 
found via molecular dynamics (MD) simulations [|12|-[l8[. 
These MD simulations do not, however, take into ac- 
count a vortex density gradient, and do not agree with 
some the varieties of plastic dynamics observed in exper- 
iments. Section IV describes the original coarse-grained 
model, first used to describe vortex avalanches H and 
vortex rivers Jl9| as seen in experiments |2Q,pl| > an d ex- 
plains the exact model used here to describe the trans- 
port experiments summarized in Section II. To eliminate 
all potentially spurious sources of noise, a completely de- 
terministic variant of the original BP model is used. To 
describe the IV experiments, we apply a shift in bound- 
ary conditions or generalized "tilt" of the vortex pile to 
represent the effect of a transport current. The resul- 
tant vortex flow represents the measured voltage. Section 
V contains the main numerical results and comparisons 
with both experiments and MD simulations. The last 
section summarizes our conclusions. 



II. BRIEF SUMMARY OF TRANSPORT 

MEASUREMENTS 

In a series of papers, Bhattacharya and Higgins 
P,E2-24| describe experiments on the nonlinear trans- 
port properties of the FLL in the anisotropic supercon- 
ductor 2H — NbSe2- In this system, the pinning is ex- 
tremely weak (jc/jo ~ 10~ 6 ), the lattice is well-formed 
and a robust "peak effect" occurs slightly below H c2 . 
Since the London length, A, is much less than the film 
thickness, the system operates in the three dimensional 
regime, where the flux line interactions decay exponen- 
tially for lengths larger than A. They use the strong mag- 
netic field dependence of the critical current (which they 
interpret in terms of a changing rigidity of the FLL) to 
explore the crossover between different type of dynamics 
including "elastic" flow, "plastic" flow, and "fluid" flow. 
In this regard, the material they study is an ideal experi- 
mental system, allowing the exploration of very different 
regimes in a well-controlled manner. 

Near the upper critical magnetic field, they observe a 
pronounced peak in the pinning force. This is referred to 
as a "peak effect" . (It should be distinguished from the 
peak in the differential resistance to be described later.) 
Equivalently, the critical current, I c , where some flux 
lines start to move, increases as the external magnetic 
field, H, increases. For I < I c the entire flux system is 
pinned, while for I > I c , some magnetic flux lines flow 
across the sample leading to the onset of a finite volt- 
age, V. We will focus mostly on results associated with 
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FIG. 1. Schematic figure illustrating the experimental IV 
measurements of Bhattacharya and Higgins in 2H — NbSe2 
at 4.2 K with varying applied field. After Fig. lb of Ref. H. 

As shown in Fig. lb of Ref. E], in the "peak regime" 
the IV characteristics of the superconductor vary enor- 
mously. For clarity, these experimental observations are 
illustrated schematically in Fig. 1. First, on increasing 
H, the critical current, J c , increases. This is the "peak 
effect" . Second, below some threshold magnetic field, the 
IV curves always rise concave upward from I c . This is 
the generic form of the IV for an FLL that is usually re- 
ported in the literature, and it is associated with an "elas- 
tic" regime. Third, when the external magnetic field en- 
ters the "peak regime" , the IV curves change drastically, 
starting as concave upwards but then bending over after 
a pronounced inflection point associated with a change 
of curvature. In the plastic regime, one obtains a charac- 
teristic S-shape IV curve. Close to onset, I c , it is concave 
upward, but then bends over as / increases further, sat- 
urating to a finite slope at large currents. Between the 
S-shaped and elastic IV curves may be a special curve 
which is always convex for I > I c . This appears at ap- 
proximately 5.8T in the experiments where the inflection 
point has moved to onset. Finally, above another mag- 
netic field value, the inflection point is at currents larger 
than those used in the experiments, and there is no satu- 
ration in the slope of the IV curves that can be observed. 
The numerical simulations presented here reproduce the 
entire progression of these curves in the plastic regime 
and their changes as a parameter in the model, repre- 
senting vortex interactions, is varied. It does not repro- 
duce the behavior in the "elastic" regime, as explained 
later. The elastic regime is also not observed in two- 
dimensional MD simulations of driven vortices near the 
onset of flow. 




0.00 



0.02 



0.04 



0.06 



0.08 



0.10 



I (A) 



FIG. 2. Schematic figures illustrating the experimental 
Differential resistance measurements corresponding to the IV 
measurements shown in Fig. 1. After Fig. Ic of Ref. [ [22| . 

Experimental measurements of the differential resis- 
tance, R — dV/dl, reveal a peak in R, corresponding 
to the inflection point for the S-shaped IV curves. The 
magnitude of this peak diminishes as the external mag- 
netic field approaches H C 2, and the position of the peak 
shifts to higher currents. Except for the largest magnetic 
fields used, the differential resistance eventually saturates 
at constant values at large enough currents, indicating 
Ohmic or fluid like behavior for sufficiently high driving. 
Also, in the plastic regime corresponding to the S-shaped 
IV curves, the usual scaling ansatz associated with dy- 
namic critical phenomena V ~ (/ — I c ) does not appear 
to hold. These experimental observations are illustrated 
schematically in Fig. 2. The numerical simulations pre- 
sented here reproduce this precise pattern of behavior for 
the differential resistance in the "peak" regime, including 
the changing position (to higher currents) and decreasing 
magnitude of the peak resistance as vortex interactions 
weaken, and the saturation to Ohmic or fluid-like behav- 
ior at high currents. 

Over a narrow range of parameters in the plastic 
regime, Bhattacharya and Higgins observed jaggedness in 
the differential resistance. This corresponds to secondary 
peaks in addition to the main peak in the differential re- 
sistance. As the external magnetic field is varied, the 
peaks can be made to appear and disappear, but for a 
given value of magnetic field and for a given sample, the 
peaks are reproducible, and thus act as "fingerprints" of 
the underlying pinning disorder. Our model also exhibits 
such IV fingerprints. 



III. SOME COMMENTS ON CONVENTIONAL 
THEORY 

The usual interpretation of the peak effect follows 
the original proposition of Pippard |[25| , that softer sys- 
tems are pinned more strongly than more rigid ones. 



This picture has been put on a formal basis by Larkin 
and Ovchinnikov's theory of collective pinning |26| . It 
was found that critical current increases as the exter- 
nal magnetic field increases near H C 2 because the FLL 
elastic moduli soften. However, collective pinning theory 
doesn't explain the shape of the IV curve outside of the 
elastic regime. Some advances have recently been made 
by Le Doussal and Giamarchi to probe the transition to 
plastic flow and other features pl| . 

MD simulations have indicated that, at least for the 
two dimensional case, the pinning forces and IV charac- 
teristics are determined by plastic deformations (tearing 
of the FLL) which fall outside the region of validity of col- 
lective pinning theory |1J,[l| . It is believed that the two 
dimensional elastic system is always unstable at onset, 
when a finite voltage first appears, and flow must take 
place in the plastic regime. These simulations can also be 
interpreted as describing two dimensional cross-sections 
of the three dimensional FLL in the plastic regime. The 
scale of the MD simulations is that of the vortex cores, 
which is much smaller than the London length, A. 

Shi and Berlinsky [0| presented a limited analytic 
treatment of the way in which "lattice defects alter the 
flow characteristics of the lattice under the influence of 
an external drive" . However both the density of lattice 
defects and the IV relation were determined by numerical 
MD simulations. These simulations showed that at high 
currents where the IV relation is linear, the defect den- 
sity drops. It increases sharply at lower currents where 
the IV relation develops an S-shaped curve. 

Further simulations have shown that the two dimen- 
sional plastic flow takes place in terms of rivers of mov- 
ing vortices separated by islands where the vortices do 
not move (£7|,g8|. This channel flow behavior has been 
observed experimentally using Lorentz microscopy [pT( . 
Previous simulations of the coarse-grained BP model re- 
vealed results consistent with these and also showed that 
the channels can form a braided river which exhibit self- 
affine (multifractal) behavior similar to fluvial braided 
rivers (see for example Ref. p9[). 

Nori and collaborators first studied, using MD simula- 
tions, flux driven into superconductors with random pin- 
ning, with the driving force solely due to the flux-density 
gradient. They elucidated many properties of the Bean 
state including the magnetic field profile, magnetization 
hysteresis loops, critical currents, vortex avalanches, and 
vortex rivers lq-|lfj|] . None of these studies using an open 
system with an overall density gradient reported the IV 
characteristics, though. 

Recently, Nori et al |16|] have simulated the IV curve as 
the FLL softens by varying the vortex interaction param- 
eter (see also Ref. |3(|). As in all previous MD studies 
of IV behavior of the FLL pa n8| , the vortices are con- 
tained in a periodic system, where they can neither enter 
nor leave. The initial condition is an ordered vortex lat- 
tice. Motion takes place via MD updates, with a uniform 
force applied to all vortices. 

The most important difference, besides the scale of the 



model, with our description of the IV experiments, is 
that, in all of these MD studies of the current-voltage 
relation, no overall vortex density gradient can develop 
owing to the periodic boundary conditions. This makes 
the vortices travel perpetually around the system and 
forces the paths to circle, which is unphysical, and does 
not occur in the real system. The actual physical sit- 
uation is an open system with flux pushed in and out, 
rather than a periodic one. The physical system adjusts 
its overall magnetic flux profile in response to applied 
forces. This is not possible in a periodic, closed system. 
Nevertheless, the basic result that the critical current in- 
creases as the vortex interactions weakens is obtained. 
However, all the IV curves measured via MD simulations 
of periodic systems, fall on top of each other, or overlap, 
at high currents, for different values of the vortex interac- 
tion strength. This is inconsistent with the experiments 
of Bhattacharya and Higgins (see Fig. 1), and reflects 
the fact that the artificially periodic system is not able 
to adjust its profile in response to applied forces. An 
even more significant, but related, difference is that the 
peak in the differential resistance grows monotonically 
and gets sharper as the resistance peak shifts to higher 
currents, the exact opposite of what happens in experi- 
ments. 

In short, the systematically changing morphology of 
the IV curves and differential resistance for a soft FLL 
has not been theoretically described. 



IV. THE MODEL 
A. Motivation 

We take a different approach to modeling flux flow in 
superconductors. In order to describe collective trans- 
port in disordered media, it is reasonable to consider 
coarse-grained models where the details of the precise 
interactions between flux lines (or vortices) are lost but 
the general effects of repulsive interaction between granu- 
lar or discrete objects, pinning, and over-damped motion 
leading to stick-slip dynamics (tearing) are preserved. 

The BP model is an interacting sand-pile model of 
vortices in a type-II superconductor, where the "sand" 
grains, representing magnetic vortices, repel each other. 
It was originally motivated by the observation of (possi- 
bly) self-organized critical [p2fl avalanches in field ramp- 
ing experiments |2(| , where the distribution of flux pack- 
ages falling into the interior coil of a hollow cylindrical 
superconductor were measured. Actually, the similarity 
between the Bean state |33| (or vortex pile) and sandpiles 
was first pointed out by de Gennes Q. Later, Vinokur, 
Feigel'man, and Geshkenbein |35) suggested that ther- 
mally induced flux creep would lead to a self-organized 
critical state in a type II superconductor, as did Tang 
fl36| . The key observation is that flux line flow always 
takes place on a flux pile which has an overall density 



gradient. This pile may be in a self-organized critical, or 
some other nontrivial state. 

In addition to describing vortex avalanches in field 
ramping experiments [g0[ 7 the BP vortex model has also 
been used to describe flux noise |37| , vortex avalanches in 
the presence of a periodic, dense array of pinning centers 
p8[ , thermally activated flux creep 39 4(J] , and magneti- 
zation loops |4l) . 



1. Coarse-graining 

Consider a transverse two-dimensional slice of a super- 
conducting slab at T = 0. The BP model B results from 
a coarse-grained description to the scale of the London 
length, A, of the microscopic dynamics, and incorporates 
the features that are essential to produce the observed 
complex behavior: Repulsive interactions between vor- 
tices, variations in the pinning potential, and variations 
in the vortex density - all at the scale of A. One can imag- 
ine imposing a grid of cells on the system; vortices in the 
model correspond to a vortex number in an extended re- 
gion of the actual physical system. Pinning in the model 
corresponds to a number of point pins in an extended cell. 
Each lattice site in the model can hold many vortices, and 
can have a different, albeit quenched, pinning potential, 
due to the underlying randomness in the positions and 
strengths of the microscopic pinning centers. The model 
allows many vortices to interact with each other while 
maintain locality (at the scale of A) in interactions. It 
is the only model of this sort that has been proposed to 
describe flux dynamics in superconductors. It enables 
numerically studies, using ordinary workstations, of the 
steady state and transient properties of systems larger 
than (500A) 2 containing tens of millions of vortices. 



B. Definition 

The BP model is defined as follows (see Fig. 3). Con- 
sider a two-dimensional honeycomb lattice |42[ where 
each cell, x, has three nearest neighbors, and is occu- 
pied by an integer number of vortices, 771(2;). The total 
energy of the vortex system includes the repulsive pair- 
wise interaction between the vortices and the attractive 
interaction of vortices with the pinning potential, V. For 
a given configuration of vortex number {777(2;)}, the total 
energy of the system is: 

H({m(x)}) = J^ JiMi)m{j) - ]T ^ pin (i)m(*) . (1) 

i,j i 

Since the model describes a system coarse-grained to the 
scale of the London length, the repulsive interactions Jy 
are short-ranged. This usually includes an on-site inter- 
action, and a weaker nearest-neighbor interaction. 

As in the microscopic case described by MD simula- 
tions, the change in the total energy of the model when 



moving a unit vortex from one site to a nearest neighbor 
site is determined. This yields the force to move a unit 
vortex from x to y, which is 



F x ^ y — Vpi n (y) Vpin 



+r[m(xl) 



(x) + [777(2;) — 771(7/) — 1] 
m(x2) — 777(7/!) — 777(7/2)] 



(2) 



As indicated in Fig. 3, the nearest neighbor cells of x 
are y, xl, and x2, and the nearest neighbors cells of y 
are x, yl, and 7/2, and < r < 1. A slightly differ- 
ent implementation of the disorder is used than before. 
The normalized pinning potential Kji n (a;) is a random 
number taken from a uniform distribution in the interval 
between zero and V max . In each time step, all cells of the 
lattice are updated in parallel. A single vortex moves 
from a cell to a neighboring cell if the force in that direc- 
tion is positive, or equivalently if the total energy of the 
system is lowered. In Eq. 2, the units of force on a vortex 
have been normalized so that the on-site term is unity. 
Thus there are two dimensionless parameters remaining, 
V max and r. 




FIG. 3. The two-dimensional honeycomb lattice. Each 
cell, x, has three nearest neighbors, and is occupied by an 
integer number of vortices, m(x). The force pushing a vor- 
tex from cell x to cell y is calculated by taking the discrete 
gradient of the sum of two potentials, one representing the re- 
pulsive interaction between vortices occupying the same and 
nearest neighbor cells, and the other representing the attrac- 
tive interaction between vortices and pinning centers. 

Many alternatives exist to handle the situation when 
more than one unstable direction appears for a vortex to 
move. In the previous implementation of the model, one 



unstable direction was chosen at random. In order to 
simplify the model and eliminate all potentially spurious 
sources of noise, the most unstable direction that has the 
largest force is chosen and the vortex moves to that site. 
This represents an extremal process. In fact the entire 
model is now completely deterministic, corresponding to 
a T = limit of the dynamics. 



C. The External Magnetic Field: Building a Vortex 

Pile 



behavior of the system. The instabilities can be elimi- 
nated by keeping track of the direction from which the 
last vortex moved onto each site and always forbidding 
a return movement backwards in that direction. Other- 
wise, backward jumps are rare and therefore disallowing 
them does not change the large scale behavior of the sys- 
tem. A similar rule applies to the boundary sites. The 
advantage of the parallel update is that it is numerically 
more efficient, than other update schemes and does not 
introduce any uncontrollable spurious effects. 



Flux lines enter the superconductor from the edges, 
pushed in by the external magnetic field. This is repre- 
sented by putting all sites on the left edge of the model 
in contact with a reservoir of vortices at some potential, 
corresponding to the external magnetic field on the left 
side of the sample. The same is done for the right side of 
the sample. For simplicity, periodic boundary conditions 
are used for the top and bottom. If the two reservoirs are 
set equal, representing embedding the sample in an exter- 
nal magnetic field, vortices enter the system generating 
the classic V-shaped flux density curve as the external 
magnetic field is increased (see below). 



1. Details about Boundary Algorithm 

In order to calculate the force on a vortex to move 
to or from a boundary site, a special algorithm must be 
used because one of the nearest neighbor sites of each 
boundary site is not on the lattice. The rule used here 
simply assumes that the "virtual" off lattice site neigh- 
boring each boundary site is occupied with an equal num- 
ber of vortices as that boundary site. More specifically, 
at the beginning of each lattice update, all of the sites 
on the boundary are set to be occupied with the same 
number of vortices. The lattice update then proceeds, 
during which vortices can move off the boundary sites 
into the system, or from the system onto the boundary 
sites, thereby changing the number of vortices occupying 
a boundary site. However, at the beginning of the next 
lattice update all of the sites on the boundary are reset to 
their original value. Through this process, vortices can 
be removed or added to the system. In general, the left 
and right boundaries are held at different values. There 
is no pinning at the boundary sites. 



2. Details about Parallel Update 

An artifact of parallel updating is the existence of local 
instabilities in which two, or more, vortices oscillate back 
and forth between neighboring sites. These local insta- 
bilities disappear if the model is coarse-grained, because 
then the neighboring sites are incorporated into a single 
one, and therefore are not important to the large scale 



D. The Transport Current: Shifting the Boundary 
Conditions 

We consider an infinite slab of finite thickness in a 
parallel applied magnetic field, carrying a current per- 
pendicular to the field. Depending on the direction of 
the applied current, the magnetic field on one side of 
the superconductor (e.g. the right hand side) will be de- 
creased and the magnetic field on the other side will be 
increased (e.g. the left hand side). Assuming the applied 
magnetic field is sufficiently strong, this corresponds just 
to changing the heights of the magnetic flux pile on the 
left and right edge, or a general kind of "tilt" of the pile. 
See, for example, Ref. E3]. Of course, the internal cur- 
rents and forces inside the superconductor will readjust 
to accommodate this new boundary condition. The flux 
lines are considered to be perfectly stiff, and described 
by a two-dimensional slice of the three dimensional slab. 
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FIG. 4. Magnetic field profile as a function of distance 
from the center of a superconducting slab in the direction 
perpendicular to both the external magnetic field and applied 
current. The long dashed line is the stable profile when an 
external field is applied, but no applied current, I. The short 
dashed, and dotted lines are the stable hysteretic profiles re- 
sulting from increasing applied current. The solid line is the 
profile just above the onset of vortex motion. 

An applied magnetic field in the geometry described 
above will produce the well known "V profile" of the 



magnetic flux density discussed, for example, in Orlando 
144J and observed in many experiments, such as those of 
Behnia, et al 45 1. This is shown in Fig. 4. The actual 



magnetic profile depends on the history of the sample 
and how magnetic field has been applied in the past. 
Applying a finite current shifts the boundary conditions 
resulting in a hysteretic profile, also shown in Fig. 4. 



1. Onset 

Eventually, as the applied current is increased further, 
a critical current is reached, where the shift of the bound- 
ary conditions is so large that steady vortex flow occurs 
down the gradient spannning the entire sample. A pro- 
file just above the critical current is also shown in Fig. 4. 
Note that this profile also retains some hysteretic proper- 
ties. For example, it has a bump corresponding to where 
the increasing and decreasing portions of the magnetic 
field profile merged as the external current was increased 
above the threshold. This bump disappears if the bound- 
ary conditions on the sample are shifted further and then 
lowered back to the previous value. 



E. Making IV Measurements 

The IV characteristic is determined by the relation be- 
tween applied current and the vortex flow, which induces 
a voltage. To our knowledge, this type of numerical mea- 
surement, made by shifting the boundary conditions on 
the Bean state, has not been investigated before. Here 
the IV characteristic is simply the relation between the 
magnitude of the shift (representing an applied current) 
and the average flow rate of vortices (representing a volt- 
age) when the critical current (tilt) is exceeded. In order 
not to confuse the reader we use the term current to refer 
to the applied electrical transport current and the term 
flow to refer to the motion of magnetic vortices. 

In general, the boundary sites in the model can be set 
to any real value, including non-integer values. This is 
because vortex number on the boundary sites describes 
the external magnetic field density at the boundary of 
the sample. However, only integer units of magnetic flux 
can enter the interior of the system, and thus only in- 
teger numbers of vortices occupy interior lattice sites. 
Obviously, the magnitude of the difference between the 
heights of the left and right boundaries, can also be set 
to non-integer values. However, if the boundary heights 
are shifted by less than whole integers, steps can appear 
in the IV data. These steps are caused by the fact that 
only whole numbers of vortices can occupy interior lat- 
tice sites. Since one vortex unit in the model represents 
many actual physical vortices, this is to some extent an 
unphysical artifact of the model, and the discreteness ef- 
fect should be less apparent in experiments. To eliminate 
this effect, all of the IV data presented in this paper were 



calculated by shifting the boundary heights by only inte- 
ger numbers of vortices. 

In simulations of the model, there is almost no back- 
ward movement of vortices. Thus, the vortex flow can 
be determined by measuring the average number of vor- 
tices moving per lattice update (the average activity). 
Furthermore, the velocity of each moving vortex is one 
lattice site per update. Thus, the experimentally mea- 
sured voltage, which is equal to the amount of moving 
flux times the velocity of that flux, is therefore propor- 
tional to the average vortex activity. 
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FIG. 5. Finite-size scaling plot of IV measurements from 
simulations of the cellular model. System sizes are shown in 
the legend. 

IV data from simulations of different size systems col- 
lapse nicely in a scaling plot, Fig. 5, if the voltage is 
measured as vortices moving per lattice update per lat- 
tice site (the average activity per site), and the current 
is measured as the average slope across the system (the 
magnitude of the height difference between the left and 
right boundaries divided by the length of the system). 
In the following sections, IV data is presented in these 
scaling units. In Fig. 5, the IV data was produced by 
repeatedly increasing the height of the left boundary by 
1 vortex and lowering the height right boundary by I 
vortex. The vortex interaction strength was r = 0.1. 

The IV data presented in the following sections was 
calculated in a similar fashion. All of the data is for sys- 
tems of size 200x400. Each IV data point was calculated 
by first shifting the left boundary height up 1 vortex and 
shifting the right boundary height down 1 vortex. As for 
Fig. 5, the lattice was updated 20000 times to eliminate 
transient behavior, and finally the lattice was updated 
another 20000 times during which the average vortex ac- 
tivity was measured. 



V. COLLECTIVE TRANSPORT BEHAVIOR OF 
THE MODEL 

The IV relation was measured for different values of 
the parameters and for different system sizes. Our main 
result is shown in Fig. 6, where the parameter V max = 5 
and the parameter r is varied. The parameter r repre- 
sents the strength of repulsion between vortices at nearest 
neighbor cells (of size A) . 
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FIG. 6. IV measurements from numerical simulations, for 
four different value of the vortex interaction strength, r. 

The first result is that as r decreases, the critical cur- 
rent, which is the slope of the pile where vortices first 
start to flow, increases. Clearly, applying an increased 
r to the pile in the steady state lowers the slope of the 
pile since formerly stable local slopes will now become 
unstable due to the increased repulsion between vortices 
at neighboring sites. Thus, fixing all other parameters, 
we can identify the parameter r in the model as a way 
of controlling the critical current, or slope of the pile. 
In the superconductor, the critical current can be con- 
trolled by the applied magnetic field. In the peak effect 
regime, it turns out that increasing the applied magnetic 
field leads to a softer FLL and thus a higher critical cur- 
rent. Therefore, the regime where the critical current 
is an increasing function of the applied magnetic field is 
represented in our model by a critical current which is 
a decreasing function of the parameter r. This is made 
evident by comparing Figs. 1,2 and Figs. 6,7. 

The second result is that, except for the largest r, all of 
the IV curves have a characteristic S shape. They start 
out at some I c increasing concave upward and then bend 
over saturating to a finite slope at large currents. The IV 
relations for a given realization of disorder, however, do 
not overlap at high currents, unlike the results obtained 
with previous MD simulations Jig ], Again, all of this 
agrees with the experimental results. 
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FIG. 7. Differential resistance measurements from numer- 
ical simulations. These results are calculated by numerical 
differentiation of the IV results shown in Fig. 6. 

In Fig. 7, the differential resistance dV/dl is calculated 
by taking numerical derivatives of the curves shown in 
Fig. 6. At the largest r value, the resistance peak is very 
large. As the parameter r decreases, the peak moves to 
higher currents and decreases in magnitude. In fact the 
same behavior is observed in experiments as the external 
magnetic field approaches H C 2- As mentioned before, the 
MD simulations give the opposite result of a peak that 
increases in magnitude as it shifts to higher currents [|16| . 

In the actual experiments, at fields below the peak ef- 
fect regime, the vortex flow is believed to be elastic and 
there is no observed peak in the differential resistance. 
As the magnetic field increases a peak in the differential 
resistance starts to develop, which then reaches a maxi- 
mum, decreasing again for larger fields. Our simulations 
appear to describe the experiments once the peak in dif- 
ferential resistance has reached its maximum. 

In order to obtain an elastic regime, we could con- 
sider a three dimensional coarse-grained model of re- 
pelling flux lines rather than point vortices representing 
a two-dimensional cross section of that system. Work is 
in progress along those lines. 

The model discussed here does not describe the be- 
havior of the superconductor in fields greater than that 
which gives the largest critical current, where the critical 
current decreases as the magnetic field increases. This 
may be due to the fact that we take the depairing cur- 
rent to be strictly infinite (i.e. j c /jo = 0) and there is no 
transition to a non-superconducting state in the model 
presented here. 

Note that the differential resistance curves for small r 
contain secondary peaks, in addition to the main peak. 
This is similar to the jaggedness or "fingerprint" found in 
experiments. This jaggedness in our results occurs in the 
filamentary channel regime, discussed below, and is due 
to filaments opening and closing as the applied current 
increases. 



A. Relation to River Morphology 

So far we have only characterized the model using mea- 
surements analogous to those that experimentalists typi- 
cally have available. However, numerical simulations can 
also provide a great deal of easily accessible information 
about the morphologies of flow patterns associated with 
different configurations. In fact, different morphologies 
of vortex flow patterns have been observed using MD 
simulations |-|lO|,[2|Jl|^,||l . To understand the na- 
ture of the differences in the shape of the IV curves, we 
have examined how the flow morphologies change as one 
increases the applied current for a system that has an S- 
shaped IV curve, and also for the singular IV curve that 
occurs at large r. 

The paths that the vortices take as they cross the sam- 
ple can be determined by measuring the average activity 
at each site. For example, Fig. 8 shows a series of grey 
scale images of the vortex flow patterns. These images 
represent "time-lapsed photographs" of the vortex activ- 
ity. The different images in Fig. 8 show the vortex flow 
patterns for the case r = 0.1 as a function of increasing 
external current, /. In these images, sites with no vor- 
tex activity are blank, sites with an activity greater than 
or equal to 0.5 (one vortex moves every other lattice up- 
date) are black, and sites with activity between and 0.5 
are indicated by grey dots with a darkness proportional 
to their activity. These patterns are fixed and do not 
change in time with fixed external driving conditions. 

The nature of the flow can be quantified by the dis- 
tribution of activity at the different lattice sites. Figure 
9 shows histograms of the activity corresponding to the 
images of Fig. 8. These histograms are constructed with 
1000 bins and normalized so that the area under the curve 
is equal to one. Note that there are peaks corresponding 
to zero activity not included in the figure. Those peaks at 
zero activity decrease in size as the current is increased. 

As can be seen in Figs. 6 and 7, the IV curve for r = 0.1 
is S-shaped. Figure 8a shows the vortex flow pattern just 
above threshold for that case. The vortex flow takes place 
only on a single filamentary string with some small side 
branches. This behavior is also evident in the histogram 
Fig. 9a, which shows a small number of isolated peaks. 
As the tilt of the pile increases the number of filaments of 
the vortex flow increases, and they begin to merge. This 
process can be seen in Figs. 8b and 8c. The corresponding 
histograms of activity shown in Figs. 9b and 9c indicate 
an increasing number of peaks, and the development of 
a continuous distribution. During this merging process, 
while the vortex flow remains filamentary, the IV curve 
remains concave upward, and the differential resistance 
rises. 

Eventually the filaments of vortex flow merge to form 
a braided river, as shown in Fig. 8d. This occurs around 
the peak in differential resistance, corresponding to a 
change in curvature of the IV relation. In this case, the 
corresponding histogram of activity (Fig. 9d) shows a 



continuous distribution of activity peaked at zero. Thus 
the peak in differential resistance signals a change in 
the underlying vortex flow morphology from filamentary 
strings to a braided river. This appears to be consistent 
with the observation of Kolton et al that at the peak in 
differential resistance "all the vortices are moving in a 
seemingly isotropic channel network with maximum in- 
terconnectivity" Jig] , 

As the tilt of the pile is increased even further, the 
vortex flow becomes spatially more and more uniform, 
as the braided river floods. This uniform flow region cor- 
responds to linear, or Ohmic behavior. These results can 
be seen in the river flow pictures in Figs. 8e and 8f. The 
corresponding histograms shows a peaked continuous dis- 
tribution, which has separated from zero activity. As the 
tilt increases further the continuous distribution of activ- 
ity obtains an increasing mean and narrowing width. The 
transition to Ohmic behavior comes about as the braided 
river floods at higher flow rates than can be supported by 
such a structure. 
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FIG. 8. Vortex flow patterns at different values of the external transport current for a vortex interaction strength of r = 0.1 
with a corresponding S-shaped IV curve. The current in each case (measured as the slope of the system) is: (a) 0.21, (b) 0.225, 
(c) 0.25, (d) 0.30, (e) 0.375, and (f) 0.45. 
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FIG. 9. Histograms of the site activity corresponding to the vortex flow patterns shown in Fig. 8. Not visible in the figures 
are the peaks at zero activity. The size of those peaks are: (a) 993, (b) 920, (c) 735, (d) 95, (e) 14, and (f) 7.6 . 



A similar progression of flow morphologies, from iso- 
lated filamentary channels, to a braided river, to a 
flooded river with uniform flow, occurs for other values 
of r with S-shaped IV curves. However, a different scene- 
rio occurs at large values of r where the IV curve is not 
S-shaped, but is instead always concave downward from 
the onset of vortex motion. The IV curve for r = 0.4 in 
Fig. 6 is an example. We will refer to this as the critical 
curve; its functional form is discussed in the next subsec- 
tion. The differential resistance for the critical curve is 
discontinuous, as can be seen in Fig. 7. 

For the critical curve, there is no filamentary region of 
vortex flow. Instead, right at the onset of vortex motion, 
the flow has the structure of a braided river. This can be 
seen in Fig. 10a, and in the corresponding activity his- 
togram in Fig. 11a which shows a continuous distribution 
of activity. (Figures 10 and 11 were produced in the same 
manner as Figs. 8 and 9, respectively.) As the tilt of the 
pile is increased even further, the flow again increases 
and becomes more spatially uniform, similar to behavior 



of the S-shaped IV case after the peak in the differential 
resistance, where the braided river floods. This behavior 
can be seen in the river flow pictures in Figs. 10b, 10c 
and lOd, and the corresponding histograms in Figs, lib, 
lie and lid. Since the critical curve has a diverging 
differential resistance at onset, this further supports our 
observation that the peak in the differential resistance is 
associated with a braided river structure. 
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FIG. 10. Vortex flow patterns at different values of the external transport current for a vortex interaction strength of r = 0.4 
with a corresponding critical IV curve. The current in each case (measured as the slope of the system) is: (a) 0.10, (b) 0.15, 
(c) 0.20, and (d) 0.40. 
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FIG. 11. Histograms of the site activity corresponding to the vortex flow patterns shown in Fig. 10. 
figures are the peaks at zero activity. The size of those peaks are: (a) 15.7, (b) 5.6, (c) 4.4, and (d) 3.3 . 
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B. Scaling of the Critical IV Curve 

The critical IV curve for r = 0.4 is well described by 

v~{i-i c y , 

with (3 = 0.6±0.1, as shown in Fig. 12. For that measure- 
ment, the critical current, I c , was measured to be 0.095, 
which was the largest current found with zero voltage. 
Since the current was sampled at values spaced by 0.005, 
the uncertainty in I c is approximately ±0.005. Varying 
the value of I c over that range changes the value of (3 
that best fits the data near the onset of vortex motion, 
and allows an estimate of the error on (3. 

The exponent [3 can be related via scaling argu- 
ments to the exponents characterizing the distribution 
of avalanches in the self-organized critical state. In the 
self-organized critical state, the average flow rate, V, of 
vortices is controlled rather than the overall slope of the 
system. The scaling argument is similar to that used in 
Corral and Paczuski [l47f to describe the transition from 



avalanche to continuous flow in the one dimensional Oslo 
model pTJ. The excess slope above the critical slope at 
onset is Am = I — I c = (AN)L 3 , where AN is the excess 
number of vortices in the system, and L is the system 
size. If vortices are added very slowly then there will be 
distinct avalanches separated by intervals of no activity. 
In that regime, superposition applies and AN ~ VL Z , 
where L z is the cutoff in the duration of the avalanches. 
In the rapidly driven regime, the avalanches overlap, and 
the excess slope becomes independent of system size, thus 
Am ~ V l 'P . These two limits can be combined into a 
single cross-over scaling function: 



AN ~ L a f{VL x ) 



(3) 



The exponent x in this expression measures the average 
duration of avalanches (t) ~ L x . Since avalanches in this 
case come about from adding an entire row of L vortices 
to the system, the average size of avalanches, measured in 
terms of the number of topplings is (s) ~ L 2 , rather than 
(s) ~ L, as in the Oslo model. Obviously, on average each 
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site in the system topples once when a row of vortices 
is added in the self-organized critical state. Using this 
result together with conservation of probability gives a 
scaling relation between x and z; 

x = z + 2 - D . 

Combining all these results gives 

x z+2- D 



a = 



D 



Using the exponent values z = 1.5 and D = 2.7 obtained 
in Ref. gives (3 = 0.6, in good agreement with the 
numerical result presented in Fig. 12. 




FIG. 12. Double logarithmic plot of voltage verses current 
minus critical current for r — 0.4. The straight line shown has 
a slope of 0.615. 



completely consistent with experimental results. How- 
ever, presumably due to the two-dimensionality of the 
model, the elastic behavior of the FLL is not reproduced. 
The success of these efforts to describe the plastic 
transport behavior of magnetic flux in superconductors 
with a coarse-grained cellular model suggests a possibly 
generic explanation for and ubiquity of plastic flow phe- 
nomena observed in superconductors. It may not depend 
on the degree of disorder or defects in an underlying mi- 
croscopic flux line lattice. Instead the behavior may be 
common to driven repulsive particle systems in a disor- 
dered media. The varieties of plastic flow behaviors in 
the model studied here result from the changing mor- 
phologies of the vortex flow pattern down a vortex den- 
sity gradient, self-organizing into distinct large scale pat- 
terns. These include isolated filaments, which merge at 
higher flow rates to give a braided river, and lead to uni- 
form flow, or Ohmic behavior, at the highest flow rates. 
The filamentary structure is associated with a concave IV 
characteristic; the braided river structure is associated 
with the peak in the differential resistance; the change to 
Ohmic behavior comes about as the braided river floods 
as it can not support a higher level of flow. 
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VI. SUMMARY 



We have studied the nonlinear current-voltage char- 
acteristics of flux flow in type II supeconductors with 
random pinning using a simple, cellular model. As in the 
physical system, vortices flow down a flux density gradi- 
ent. Because the coarse-grained model does not include 
any degrees of freedom at scales much smaller than the 
London length, A, it can not describe structural ordering 
or disordering of the flux line lattice. Despite this fact, 
simulations of the model reproduce many of the empir- 
ically observed transport characteristics of plastic flux 
flow in type II superconductors. 

In particular, our results reproduce many of the fea- 
tures attributable to plastic flow of the FLL in the "peak 
regime" . By weakening the vortex interaction strength, 
we find an increase of the critical current and a falling of 
the magnitude of the peak in the differential resistance. 
The model also exhibits IV fingerprints, and crossover 
to Ohmic or linear behavior at high currents. Also, the 
IV curves for different vortex interaction strengths do not 
merge at large external currents. All of these features are 
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